Now we load the packages we will require along the analysis.

source("sdrtools.R")
## Loading required package: distrom
## Loading required package: Matrix
## Loading required package: gamlr
## Loading required package: parallel
## Loading required package: dirmult
## Loading required package: MASS
## Loading required package: glmnet
## Loading required package: foreach
## Loaded glmnet 2.0-16
## Loading required package: pracma
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:Matrix':
## 
##     expm, lu, tril, triu

First we load the data and check for possible variables without expression. Data comes in three different files. One with the compositional data (relative abundance), one with the total number of reads for each sample and one with the labels.

compositions = read.csv('./data/HMP/HMPdataL2.txt',header=FALSE)
totals = read.table("./data/HMP/HMPdataL2reads.txt")
labels = read.table("./data/HMP/HMPlabels.txt")
HMPdata = checkData(compositions,totals,labels)
X.counts = HMPdata$counts
Y = HMPdata$labels

sites = c('Stool','Saliva','Skin','Nasal','Vagina')
groups = NULL
for (n in 1:length(Y)){
  groups[n] = sites[Y[n]]
}
groups = as.factor(groups)

We start with visualization. We first analize the data applying widely used tools derived from pairwise distance matrices. In particular, we use PCoA as implemented in APE and MDS as implemented in VEGAN.

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-5
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:pracma':
## 
##     procrustes
library(ape)
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
symbols = c("triangle-up-dot","circle-open","diamond-open","cross","square-open")

logX = log(X.counts+1)
Dmtx = vegdist(logX)
res = pcoa(Dmtx)

pcoa2plot = data.frame(Y=groups,X1=res$vectors[,1],X2=res$vectors[,2])
ggplot(pcoa2plot, aes(x=X1, y=X2, color=groups, shape=groups)) + 
  geom_point(size=3, alpha=0.8) + xlab('PCoA-1') + ylab('PCoA-2')

nmds = metaMDS(X.counts,autotransform = FALSE,trace=0)

nmds2plot = data.frame(Y=groups,X1=nmds$points[,1],X2=nmds$points[,2])
ggplot(nmds2plot, aes(x=X1, y=X2, color=groups, shape=groups)) + 
  geom_point(size=3, alpha=0.8) + xlab('MDS-1') + ylab('MDS-2')

Let’s take a look first to the full-rank reduction:

normal.fit = sdr4normal.fit(X.counts,Y,4)
projLN = sdr4normal.project(normal.fit,X.counts)

par(mfrow=c(2,2))
boxplot(projLN[,1]~groups, main="Projection onto Normal-SDR-1",ylab = "SDR-1")
boxplot(projLN[,2]~groups, main="Projection onto Normal-SDR-2",ylab = "SDR-2")
boxplot(projLN[,3]~groups, main="Projection onto Normal-SDR-3",ylab = "SDR-3")
boxplot(projLN[,4]~groups, main="Projection onto Normal-SDR-4",ylab = "SDR-4")

Now we use SDR with the normal model on the log-contrast transformation. We first estimate the best dimension for the reduction:

d = testDim4sdr(X.counts,Y,"normal",alpha=0.05)
print(d)
## [1] 2

We now compute the estimate corresponding to the estimated dimension:

normal.fit = sdr4normal.fit(X.counts,Y,d,lambda=1.0)
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Lambda: 1  nr.var: 44
projLN = sdr4normal.project(normal.fit,X.counts)

resuLN = data.frame(Y=groups,X1=projLN[,1],X2=projLN[,2])
ggplot(resuLN, aes(x=X1, y=X2, color=groups, shape=groups)) + 
  geom_point(size=3, alpha=0.8) + xlab('N-SDR-1') + ylab('N-SDR-2')

Let’s repeat the analysis using the multinomial model. Let’s take a look first to the full-rank reduction:

mn.fit = sdr4multinomial.fit(X.counts,Y,4)
projMN = sdr4multinomial.project(mn.fit,X.counts)

par(mfrow=c(2,2))
boxplot(projMN[,1]~groups, main="Projection onto Multinomial-SDR-1",ylab = "SDR-1")
boxplot(projMN[,2]~groups, main="Projection onto Multinomial-SDR-2",ylab = "SDR-2")
boxplot(projMN[,3]~groups, main="Projection onto Multinomial-SDR-3",ylab = "SDR-3")
boxplot(projMN[,4]~groups, main="Projection onto Multinomial-SDR-4",ylab = "SDR-4")

Now we test for the dimension and repeat the analysis, including variable selection:

d = testDim4sdr(X.counts,Y,"multinomial")
print("Estimated dimension is:")
## [1] "Estimated dimension is:"
d
## [1] 2
mn.fit = sdr4multinomial.fit(X.counts,Y,d,lambda=1.5)
## [1] 0
## [1] 0
projMN = sdr4multinomial.project(mn.fit,X.counts)

resuMN = data.frame(Y=groups,X1=projMN[,1],X2=projMN[,2])
ggplot(resuMN, aes(x=X1, y=X2, color=groups, shape=groups)) + 
  geom_point(size=3, alpha=0.8) + xlab('MN-SDR-1') + ylab('MN-SDR-2')

Now we repeat the analysis using a PGM model. We start we the full-rank estimate before testing for the dimension:

pgm.fit = sdr4pgmR.fit(X.counts,Y,4)
projPGM = sdr4pgmR.project(pgm.fit,X.counts)

par(mfrow=c(2,2))
boxplot(projPGM[,1]~groups, main="Projection onto PGM-SDR-1",ylab = "SDR-1")
boxplot(projPGM[,2]~groups, main="Projection onto PGM-SDR-2",ylab = "SDR-2")
boxplot(projPGM[,3]~groups, main="Projection onto PGM-SDR-3",ylab = "SDR-3")
boxplot(projPGM[,4]~groups, main="Projection onto PGM-SDR-4",ylab = "SDR-4")

We now test for the dimension of the reduction:

d = testDim4sdr(X.counts,Y,"PGM")
print(d)
## [1] 3
pgm.fit = sdr4pgmR.fit(X.counts,Y,d,2.5)
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Lambda: 2.5  nr.var: 32
projPGM = sdr4pgmR.project(pgm.fit,X.counts)

resuPGM = data.frame(Y=groups,X1=projPGM[,1],X2=projPGM[,2],X3=projPGM[,3])
library(plotly)
myscene = scene4plot3d("FPGM")
plot_ly(resuPGM,x=resuPGM$X1,y=resuPGM$X2,z=resuPGM$X3,type="scatter3d",mode="markers", marker = list(size=3,alpha=0.8), symbols = symbols, color = resuPGM$Y,symbol = resuPGM$Y) %>%
layout(scene=myscene)

To test for overall association, we fit a linear model with the reduction and perform a MANOVA analysis:

library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:pracma':
## 
##     logit
mod = lm(pgm.fit$proj~groups)
Manova(mod)
## 
## Type II MANOVA Tests: Pillai test statistic
##        Df test stat approx F num Df den Df    Pr(>F)    
## groups  4    2.1351   187.48     16   2620 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We have already used variable selection when fitting the reductions prior to visualization. Let us point out that the variable selection induced here is different form the one in Taddy (Annals of Appled Statistics, 2015), Distributed Multinomial Regression. In particular, our algorithms induce structured variable selection, which allows for effective selection of whole rows of the coefficient matrix. Let’s see this in the example with the HMP data set:

mn.fit = dmr(cl=NULL,covars=get_fyZ(Y),counts=X.counts)
(bhat_taddy = t(as.matrix(coef(mn.fit)))[,-1])
##             1          2          3          4
## 1   2.2092849  1.5736738  2.9216353  0.6557727
## 2   0.0000000  0.0000000  6.4297195  4.1292305
## 3  -1.8869328  0.0000000  2.5322352  2.7279557
## 4   0.0000000  0.0000000  6.4149883  1.7997768
## 5   3.2567244  2.2989752  0.5560673 -0.3741656
## 6   0.0000000  0.0000000  6.4161635  2.0288958
## 7   0.0000000  0.8267888  3.4253003  2.7585980
## 8   0.0000000  0.6745854  4.2610555  2.1534522
## 9  -1.1007975 -0.7938358 -1.4841515 -1.2113223
## 10 -0.6250184  3.5545352  2.2113717  0.0000000
## 11 -0.3936362  4.7154943  4.1747611  0.0000000
## 12  0.0000000  0.0000000  7.3767519  0.0000000
## 13  4.5442736  0.0000000  0.0000000  0.0000000
## 14  0.0000000  0.0000000  7.3730903  0.0000000
## 15  0.9086664  3.0347895  2.7812210  1.9923867
## 16  0.0000000  5.7515207  0.0000000  0.0000000
## 17 -0.2287752  4.6095047  1.2733539  0.0000000
## 18  0.0000000  2.9456868  0.0000000  0.0000000
## 19  0.0000000  5.3991104  0.0000000  0.0000000
## 20  1.7871524  0.7100284  0.2335475 -1.9988618
## 21 -1.5246576  0.7812930  4.6838665  3.4680050
## 22  5.0784664  0.0000000  0.0000000  0.0000000
## 23  0.0000000  0.0000000  7.3700755  0.0000000

It is clear that Taddy’s estimator does not induce any variable selection but only sparsity of the coefficient matrix:

print("The selected variables are:")
## [1] "The selected variables are:"
(idxsel_taddy = which(apply(bhat_taddy,1,Norm)>0.0000001))
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

Now, applying our method for the multinomial we get:

mn.fit = sdr4multinomial.fit(X.counts,Y,4,lambda=1.5)
## [1] 0
## [1] 0
(bhat = mn.fit$bhat)
##             [,1]        [,2]       [,3]       [,4]
##  [1,]  0.0000000  0.00000000  0.0000000  0.0000000
##  [2,]  0.0000000  0.00000000  0.0000000  0.0000000
##  [3,] -2.9929572  0.53502482 -1.7842172  2.1175293
##  [4,]  0.0000000  0.00000000  0.0000000  0.0000000
##  [5,]  2.5367184 -0.64974340  0.4076681 -0.9738687
##  [6,]  0.0000000  0.00000000  0.0000000  0.0000000
##  [7,] -0.5904832  1.57979514 -1.0501777  1.2966072
##  [8,]  0.0000000  0.00000000  0.0000000  0.0000000
##  [9,] -1.9061826 -3.32734946 -2.6867410 -1.7174680
## [10,] -0.8959869  0.48030117  1.5388183 -0.1782029
## [11,]  0.0000000  0.00000000  0.0000000  0.0000000
## [12,]  0.0000000  0.00000000  0.0000000  0.0000000
## [13,]  0.0000000  0.00000000  0.0000000  0.0000000
## [14,]  0.0000000  0.00000000  0.0000000  0.0000000
## [15,]  0.3171647  1.20450541  1.3153471  1.6055153
## [16,]  0.0000000  0.00000000  0.0000000  0.0000000
## [17,] -0.6942391  0.23583900  1.8724128 -0.2832400
## [18,]  0.0000000  0.00000000  0.0000000  0.0000000
## [19,] -0.7629803  0.32464703  1.9785345 -0.3822164
## [20,]  1.5663081 -0.28574626 -0.5225408 -0.6549143
## [21,]  0.0000000  0.00000000  0.0000000  0.0000000
## [22,]  3.4226382 -0.09727345 -1.0691042 -0.8297413
## [23,]  0.0000000  0.00000000  0.0000000  0.0000000
print("The selected variables are:")
## [1] "The selected variables are:"
(idxsel = which(apply(bhat,1,Norm)>0.0000001))
##  [1]  3  5  7  9 10 15 17 19 20 22

We now discuss performance in prediction. For this we run a 10-fold cross-validation experiment using standard sampling to generate the CV partitions. (To reduce computing time, we have set a fixed value of the regularization parameter for all the aprtitions. Thus, the resulting error rates can be slightly higher than the minimum attainable)

set.seed(100)
kfold = 10
parts = sample(1:kfold,length(Y),replace=TRUE)

Xtrain = list()
Xtest = list()
Ytrain=list()
Ytest=list()

for (k in 1:kfold){
  Xtrain[[k]] = X.counts[parts!=k,]
  Xtest[[k]] = X.counts[parts==k,]
  Ytrain[[k]] = Y[parts!=k]
  Ytest[[k]] = Y[parts==k]
}

errores.mn = numeric(kfold)
errores.ln = numeric(kfold)
errores.p = numeric(kfold)
errores.pgm = numeric(kfold)
capture.output({
for (k in 1:kfold){
  print(k)
  # reduce using sdr-normal
  fit.ln = sdr4normal.fit(Xtrain[[k]],Ytrain[[k]],dim=2,lambda = 15.0)
  classfit = glmnet(fit.ln$proj,Ytrain[[k]],family="multinomial")
  proj.ln = sdr4normal.project(fit.ln,Xtest[[k]])
  probs = predict(classfit,proj.ln,s=0.0,response="class")
  yhat = apply(probs,1,which.max)
  errores.ln[k] = mean(yhat!=Ytest[[k]])  
  
  # reduce using sdr-multinomial
  fit.mn = sdr4multinomial.fit(x=Xtrain[[k]],y=Ytrain[[k]],lambda=1.5)
  classfit = glmnet(fit.mn$proj,Ytrain[[k]],family="multinomial")
  proj.mn = sdr4multinomial.project(fit.mn,Xtest[[k]])
  probs = predict(classfit,proj.mn,s=0.0,response="class")
  yhat = apply(probs,1,which.max)
  errores.mn[k] = mean(yhat!=Ytest[[k]])
  
  # reduce using sdr-PGM
  fit.pgm = sdr4pgmR.fit(x=Xtrain[[k]],y=Ytrain[[k]])
  classfit = glmnet(fit.pgm$proj,Ytrain[[k]],family="multinomial")
  proj.pgm = sdr4pgmR.project(fit.pgm,Xtest[[k]])
  probs = predict(classfit,proj.pgm,s=0.0,response="class")
  yhat = apply(probs,1,which.max)
  errores.pgm[k] = mean(yhat!=Ytest[[k]])
  
  
  # benchmark without reduction
  auxfit = cv.glmnet(Xtrain[[k]],Ytrain[[k]],family="multinomial")
  classfit = glmnet(Xtrain[[k]],Ytrain[[k]],family="multinomial")
  probs = predict(classfit,Xtest[[k]],s=auxfit$lambda.1se,response="class")
  yhat = apply(probs,1,which.max)
  errores.p[k] = mean(yhat!=Ytest[[k]])  
}
})
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -91); Convergence for 91th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -83); Convergence for 83th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -96); Convergence for 96th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -88); Convergence for 88th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -84); Convergence for 84th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -99); Convergence for 99th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -97); Convergence for 97th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -98); Convergence for 98th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -95); Convergence for 95th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -96); Convergence for 96th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -99); Convergence for 99th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -83); Convergence for 83th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -93); Convergence for 93th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -98); Convergence for 98th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -96); Convergence for 96th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -93); Convergence for 93th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -94); Convergence for 94th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -96); Convergence for 96th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -92); Convergence for 92th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -93); Convergence for 93th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -82); Convergence for 82th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned

## Warning: from glmnet Fortran code (error code -82); Convergence for 82th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -77); Convergence for 77th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -84); Convergence for 84th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -83); Convergence for 83th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -84); Convergence for 84th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -86); Convergence for 86th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -92); Convergence for 92th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -60); Convergence for 60th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -77); Convergence for 77th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in cond/sqrt(cond.norm2): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## Warning: from glmnet Fortran code (error code -93); Convergence for 93th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -89); Convergence for 89th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
## Warning: from glmnet Fortran code (error code -90); Convergence for 90th
## lambda value not reached after maxit=100000 iterations; solutions for
## larger lambdas returned
##  [1] "[1] 1"                  "Lambda: 15  nr.var: 6 "
##  [3] "[1] 0"                  "[1] 0"                 
##  [5] "[1] 2"                  "Lambda: 15  nr.var: 6 "
##  [7] "[1] 0"                  "[1] 0"                 
##  [9] "[1] 3"                  "Lambda: 15  nr.var: 6 "
## [11] "[1] 0"                  "[1] 0"                 
## [13] "[1] 4"                  "Lambda: 15  nr.var: 6 "
## [15] "[1] 0"                  "[1] 0"                 
## [17] "[1] 5"                  "Lambda: 15  nr.var: 6 "
## [19] "[1] 0"                  "[1] 0"                 
## [21] "[1] 6"                  "Lambda: 15  nr.var: 6 "
## [23] "[1] 0"                  "[1] 0"                 
## [25] "[1] 7"                  "Lambda: 15  nr.var: 6 "
## [27] "[1] 0"                  "[1] 0"                 
## [29] "[1] 8"                  "Lambda: 15  nr.var: 6 "
## [31] "[1] 0"                  "[1] 0"                 
## [33] "[1] 9"                  "Lambda: 15  nr.var: 6 "
## [35] "[1] 0"                  "[1] 0"                 
## [37] "[1] 10"                 "Lambda: 15  nr.var: 6 "
## [39] "[1] 0"                  "[1] 0"
(apply(cbind(errores.p,errores.ln,errores.mn,errores.pgm),2,median)) 
##   errores.p  errores.ln  errores.mn errores.pgm 
##  0.13049645  0.10633803  0.10319149  0.07968631